 function SAT_SST_Step_03_get_scaling(P)
        
    Out      = read_data(P);
    
    % *********************************************************************
    % Compute output arguments
    % *********************************************************************
    for ct_x = 1:36
        for ct_y = 1:18
            
            Out_grid = prepare_grid_data(Out,ct_x,ct_y,P);
            
            if ~isempty(Out_grid)
                
                clear('T_a','dT_a_dt')
                T_a     = BB98d_forcing_month2high_reso(Out_grid.AT_itp);
                T_a     = T_a - nanmean(T_a(1:100));
                dT_a_dt = [(T_a(2:end) - T_a(1:end-1)) ./ (86400/2); 0];

                % Fit parameters assuming no seasonality ------------------
                % x_fit = BB98d_fit_parameters(T_a,dT_a_dt,T_o,l_use_loss,do_season,do_spec)
                x_fit = BB98d_fit_parameters(T_a,dT_a_dt,Out_grid.SST,Out_grid.l_use_loss,P);
                P_fit = BB98d_predict(T_a,dT_a_dt,x_fit);
            
                Out_grid.AT_cowtan = Out_grid.AT * (0.86 - 0.25*Out_grid.lf/100);
                Out_grid.AT_scl    = P_fit.T_o_m;
                Out_grid.AT_scl(isnan(Out_grid.AT)) = nan;

                if P.do_season == 0
                    Parameters(ct_x,ct_y,1:3) = x_fit;
                elseif P.do_season == 1
                    Parameters(ct_x,ct_y,1:9) = x_fit;
                end
                Save_data{ct_x,ct_y} = Out_grid;
            else
                if P.do_season == 0
                    Parameters(ct_x,ct_y,1:3) = nan;
                elseif P.do_season == 1
                    Parameters(ct_x,ct_y,1:9) = nan;
                end
                Save_data{ct_x,ct_y} = [];
            end
            
        end
    end
    
    % *********************************************************************
    % Save scaling factors
    % *********************************************************************
    file_save = [SAT_SST_func_IO('data',P),'/BB98d_parameters_starting_year_',num2str(P.yr_learn(1)),'_',P.date,'.mat'];
    save(file_save,'Parameters','Save_data','-v7.3');
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function Out_grid = prepare_grid_data(Out,ct_x,ct_y,P)

    l = Out.lon > (ct_x*10-10)  & Out.lon <= (ct_x*10) & ...
        Out.lat > (ct_y*10-100) & Out.lat <= (ct_y*10-90);

    if nnz(l) > 0
        disp([ct_x ct_y])
        Out_grid = pre_process_data(Out,l,P);
    else
        Out_grid = [];
    end
end

function Out_grid = pre_process_data(Out,l,P)

    Out_grid.N   = nnz(l);

    Out_grid.lf  = nanmean(Out.land_fraction(l));

    AT  = CDC_pnt2glbmean(Out.raw_SAT_anm(l,:,:),Out.lon(l),Out.lat(l));
    SST = CDC_pnt2glbmean(Out.raw_SST_anm(l,:,:),Out.lon(l),Out.lat(l));

    yrs = P.yr_learn - 1849;
    AT  = AT(:,yrs);  AT  = AT(:);
    SST = SST(:,yrs); SST = SST(:);

    if nnz(~isnan(AT + SST)) > 240

        % Subset years where both AT and SST has values ---------------
        temp = reshape(AT+SST,12,numel(P.yr_learn));
        l_yr = any(~isnan(temp),1);
        l_yr(find(l_yr,1):find(l_yr,1,'last')) = 1;
        yrs  = P.yr_learn;
        Out_grid.yrs  = yrs(l_yr);

        % Linear interpolate if missing values ------------------------
        l_yr = repmat(l_yr,12,1);
        Out_grid.AT   = AT(l_yr);
        Out_grid.SST  = SST(l_yr);

        Out_grid.AT_itp  = CDC_interp1(Out_grid.AT,[],1);
        Out_grid.SST_itp = CDC_interp1(Out_grid.SST,isnan(Out_grid.AT),1);

        Out_grid.l_use_loss  = ~isnan(Out_grid.AT) & ~isnan(Out_grid.SST);
    else
        Out_grid = [];
    end
end


function Out = read_data(P)

    % *********************************************************************
    % Load data
    % *********************************************************************
    file_load = [SAT_SST_func_IO('data',P),'/Raw_SAT_SST_',P.date,'.mat'];
    load(file_load,'raw_SST_anm','lon','lat','land_fraction');

    file_load = [SAT_SST_func_IO('data',P),'/Raw_SAT_anm_',P.date,'.mat'];
    load(file_load,'raw_SAT_anm')

    % *********************************************************************
    % Subset long stations
    % *********************************************************************
    % GLB_mod_subset_long_stations;

    raw_SST_anm     = raw_SST_anm - nanmean(raw_SST_anm(:,:,[1982:2014]-1849),3);
    raw_SAT_anm     = raw_SAT_anm - nanmean(raw_SAT_anm(:,:,[1982:2014]-1849),3);

    Out.raw_SST_anm = raw_SST_anm;
    Out.raw_SAT_anm = raw_SAT_anm;
    Out.lon         = lon;
    Out.lat         = lat;
    Out.land_fraction = land_fraction;
    
end
